On the variability of the Bering Sea Cold Pool and implications for the biophysical environment

The Bering Sea experiences a seasonal sea ice cover, which is important to the biophysical environment found there. A pool of cold bottom water (<2°C) is formed on the shelf each winter as a result of cooling and vertical mixing due to brine rejection during the predominately local sea ice growth. The extent and distribution of this Cold Pool (CP) is largely controlled by the winter extent of sea ice in the Bering Sea, which can vary considerably and recently has been much lower than average. The cold bottom water of the CP is important for food security because it delineates the boundary between arctic and subarctic demersal fish species. A northward retreat of the CP will likely be associated with migration of subarctic species toward the Chukchi Sea. We use the fully-coupled Regional Arctic System Model (RASM) to examine variability of the extent and distribution of the CP and its relation to change in the sea ice cover in the Bering Sea during the period 1980–2018. RASM results confirm the direct correlation between the extent of sea ice and the CP and show a smaller CP as a consequence of realistically simulated recent declines of the sea ice cover in the Bering Sea. In fact, the area of the CP was found to be only 31% of the long-term mean in July of 2018. In addition, we also find that a low ice year is followed by a later diatom bloom, while a heavy ice year is followed by an early diatom bloom. Finally, the RASM probabilistic intra-annual forecast capability is reviewed, based on 31-member ensembles for 2019–2021, for its potential use for prediction of the winter sea ice cover and the subsequent summer CP area in the Bering Sea.


Introduction
The Cold Pool (CP) is the region of the Bering Sea shelf where bottom water is < 2˚C throughout the summer [1]. Cooling and seasonal sea ice formation in winter results in the formation of this cold, salty and dense water mass [1]. The freezing of surface water across the northern  [18]) and examine the linkage between CP area and sea ice area. The distribution of sea ice is shown to exert a strong control on the size of the CP and it is also found to have an important effect on the timing of phytoplankton blooms on the shelf. We also analyze results from the RASM ensemble forecasts of sea ice and bottom water temperature in the Bering Sea, for their potential usefulness to stakeholders in the region.

Model description, methods, and validation
RASM has been developed to address some of the limitations of global Earth System models, including the representation of Arctic relevant processes within and coupling among model components and high spatio-temporal resolution. It is a fully-coupled regional Earth system model with components including, atmosphere, ocean, sea ice, marine biogeochemistry, land hydrology and a river routing scheme. All the components are coupled using the flux coupler of Craig et al. [19]. The RASM domain includes all the sea-ice covered ocean in the Northern Hemisphere, the Arctic river drainage, and large-scale atmospheric weather patterns. This domain includes all oceanic areas of the northern hemisphere from 90˚N to 55˚N and most of the North Pacific down to 30˚N. We use it here as an optimal tool to investigate the Cold Pool and related bio-physical processes in detail.
The RASM components and their configurations are described in Table 1. The ocean biogeochemistry (BGC) component in RASM is a medium-complexity Nutrients-Phytoplankton-Zooplankton-Detritus (NPZD; [20]) model. The sea ice model uses a BGC parameterization, in which biological activity is distributed throughout the ice column [21]. RASM has been in development by a dedicated team for over a decade. Each component of this system model, as well as the fully-coupled system, has been validated and simulation results evaluated extensively [22][23][24][25][26][27][28][29][30].
The RASM historical  simulation results analyzed here were produced after a 78-year spinup, which started with no sea-ice and the Polar Science Center Hydrographic Climatology (PHC) 3.0 [31] climatological ocean temperature and salinity at rest and was forced with the Coordinated Ocean-sea ice Reference Experiments version 2 (CORE2) reanalysis [32]. The RASM 6-month 31-member ensemble forecast results presented here were initialized on February 1 of each forecast year. The forcing along the lateral boundaries and in the upper atmosphere, as well as atmospheric initial conditions, were derived from the National Centers for Environmental Predictions (NCEP) global Coupled Forecast System (CFS) version 2 (CFSv2) operational 9-month forecasts initialized at 0000Z each day of the preceding month. The ocean and sea ice initial conditions for these forecasts have been produced from a separate fully-coupled RASM hindcast simulation (without BGC components), which started in September 1979 and has been updated through January 2021 using surface atmospheric forcing from the NCEP global CFS Reanalysis (CFSR). The sea ice representation by RASM has been recently evaluated and validated by Watts et al. [33]. This work utilized spatial analysis metrics, such as the integrated ice edge error, Brier score, and spatial probability score to examine the historical representation  of sea ice extent, volume, and thickness. Results from RASM, as well as 12 other models which participated in phase 6 of the Coupled Model Intercomparison Project (CMIP6), were compared with the available satellite observations of sea ice parameters. Of the 13 models that were examined, RASM had the best spatial probability score and lowest integrated ice edge error for representation of the Bering Sea ice extent during the months of June to November (see Table 3  Results from the predecessor to this present configuration of RASM BGC was reported on in Jin et al. [28]. This work found that moving to a higher spatial resolution and adding an improved mixed layer depth scheme improved simulation of nitrate concentrations and primary production, especially in areas with sharp bathymetric gradients such as along shelf breaks in the Western Arctic. Their 9-km model configuration (R9km), which produced results with the lowest biases in the Jin et al. [28] study, is the predecessor to the updated model version used and analyzed here. Some of the key updates we have made include upgrading CICE to version 6, as well as utilizing the fully-coupled configuration of RASM.
More recently, results from the current BGC component of RASM have been examined by Frants et al. [34] and Clement Kinney et al. [24]. In the Frants et al. [34] study, it was found that RASM realistically represents the phytoplankton bloom of pelagic diatoms found under sea ice in a similar place and time to observations in the Chukchi Sea from 2011. However, the timing of the modeled bloom was approximately 2 weeks earlier, due to an earlier retreat of the sea ice over the northern Chukchi Sea, likely caused by discrepancies in the predicted wind field over the region. In addition, the modeled values of chl-a underestimated the in situ measurements, especially in terms of the maximum chl-a concentrations, likely due to the underrepresentation of the mesoscale physical processes, e.g. eddies or vertical mixing along the ice edge, and associated nutrient concentrations potentially captured by point measurements. The recent Clement Kinney et al. [24] study was motivated by limited observations of phytoplankton blooms beneath Arctic sea ice and the lack Arctic-wide primary production observational estimates beneath the ice. Using the fully-coupled RASM with marine biogeochemistry it was estimated that 63%/41% of the total primary production in the Arctic occurs in waters covered by sea ice that is �50%/�85% concentration. This means that the majority of total primary production in the Arctic may be missing from current observational estimates based on satellite data. The total primary production in the Arctic was found to be increasing at a rate of 5.2% per decade during 1980-2018. Increased light transmission, due to the removal of sea ice, more extensive melt ponds, and thinner sea ice, was implicated as the main cause of increasing trends in primary production.

Results
The Bering Sea bathymetry, relevant place names, and area for CP calculation are shown in March (691,000 km 2 ) and a minimum in October (180,000 km 2 ). Because recent CP observational work [7] focuses on the summer bottom water temperatures, we will also focus on the July CP values. The annual cycle of individual years is shown as black lines with a large range

PLOS ONE
On the variability of the Bering Sea Cold Pool and implications for the biophysical environment of variability (1 standard deviation during July is +/-118,000 km 2 ). Two years exhibit July CP area values outside of 2 standard deviations; these are 2012 (649,000 km 2 or 290,000 km 2 above the mean) and 2018 (113,000 km 2 or 246,000 km 2 below the mean).     [7] RASM indicates that, during the period of satellite data, 2012 was the year of maximum sea ice area (713,000 km 2 ) and 2018 was the minimum (139,000 km 2 ). RASM also shows another minima during 2001 that is also less than 2 s.d. below the mean. This minima was associated with the third lowest July CP area (Fig 4) and is discussed further in the Discussion.
Monthly mean values of ice area and CP area (over the years 1980-2018) show a significant (p-value < 0.01) correlation with a value of 0.79. Even after removing the annual cycle, the correlation is still significant with a value of 0.60.  Fig 9B and 9C), RASM is representing the spatial distribution very well, but the extent does not reach as far south as the data from NSIDC.
Modeled primary production is strongly affected by the presence or absence of sea ice and the timing of its formation and melt. Fig 10 shows   bloom with values in excess of 100 mg C / m 3 / day, whereas June 2012 shows a dwindling bloom. This can be contrasted with 2018 where May has weak values on the central and southeastern shelf, but considerably higher production during June. Vertical profiles at 60.5˚N, 175˚W (location is shown as small black circles in Fig 10) of temperature, salinity, and density from RASM reveal differences in the stratification of the water column during these two years (Fig 11). Beginning with May, there is stratification of the water column in 2012 due to salinity because of sea ice melt; however, during 2018 there was very little sea ice formation/melt and therefore no stratification of the water column. Moving on to June, we see the development of thermal stratification during both 2012 and 2018. The thermal stratification in June 2018 allows for the development of the phytoplankton bloom seen in Fig 10. This June bloom in 2018 is later than the peak of the bloom in 2012, which occurred in May.

Discussion
The extreme reduction in sea ice formation during 2018 resulted in the smallest CP area over the model timeseries from 1980-2018 (Fig 4).   The low ice area is largely due to the fact that the Bering Sea experienced strong winds out of the south in winter 2018 that restricted the typical southward expansion of sea ice toward the shelf break [7]. It is interesting to note that the winter sea ice minimum was reached only 6 years after the timeseries maximum in 2012, which reflects the strong dependence of the winter sea ice cover on the wind field in the Bering Sea.
We have focused much attention on the dramatically low CP and sea ice area during 2018.  (Fig 4) was not preceded by such anomalously low sea ice area as in 2001. However, the sea ice melted much earlier than the mean in April and May 2003 (see purple line in Fig 8), which may have contributed to the reduced CP area later in the summer.
Interannual variability in sea ice presence affects not only the formation of the CP, but also the timing of spring phytoplankton blooms on the shelf. Variability in primary production appears to follow the Oscillating Control Hypothesis originally developed by Hunt et al. (2002) for the southeastern Bering Sea. The hypothesis states that early ice retreat will lead to a late bloom, while late ice retreat leads to an early bloom. During a heavy ice year, like 2012, the

PLOS ONE
bloom happens immediately after sea ice melt in May when low surface salinity helped stratify the water column. On the other hand, during a light ice year, like 2018, the bloom across the central shelf did not occur until June after thermal stratification had developed. This has implications for the pelagic versus benthic ecosystems. In the current regime, where blooms occur early in the season, much of the production sinks to the bottom due to limited grazing by zooplankton. The sinking of this production helps to support a rich benthic community in the Bering Sea, as well as the rest of the Pacific Arctic. Food availability to the benthos is one of the largest drivers for both the benthic community as a whole [38,39] and to specific species such as Macoma calcarea, a common bivalve in the area [40]. Bottom water temperature, in this case the presence of the CP, has also been shown as a physical driver for some benthic organism populations [40]. Therefore, years with extensive sea ice are expected to likely favor the benthos since much of the production will continue to sink to the bottom without being consumed and grazed upon by zooplankton. However, in years with little sea ice and warmer temperatures there could be a shift to a more pelagic community, as zooplankton growth, abundance, and grazing could increase [41], and zooplankton species that are advected into the system may be able to survive better in these warmer waters [42], thus weakening the notable pelagic-benthic coupling of the region.
Due to the high degree of interannual variability in sea ice area and CP area on the Bering Sea shelf, it would be useful to have forecasts of these fields skillful enough to provide stakeholders with an estimate of conditions for the upcoming months. Such an objective has motivated the development of RASM intra-annual (i.e., up to 6 months) ensemble forecasts. Every month since January 2019, RASM has been used to generate (28)(29)(30)(31); depending on the number of days in the month) probabilistic forecasts of the sea ice conditions into the future, by dynamically downscaling the global CFSv2 operational 9-month forecasts (https://nps.edu/ web/rasm/welcome). While these products are not meant for commercial or operational use, they provide useful and often more realistic insights on conditions at a regional to local scale, in addition to what is available from global forecasts. Fig 12A shows  Similarly, Fig 13A shows the prediction of sea ice extent for April 15, 2020, initialized on February 1, 2020. As in Fig 12A, a very large spread exists between the minimum and maximum ice extents predicted from the ensemble of 31 members. In both years the minimum ensemble member shows an ice edge north of St. Lawrence Island, while the maximum ensemble member shows ice approximating the shelf break in the western and central Bering Sea. Fig  13B shows the probability of sea ice area >15% on April 15, 2020 based on a 31-member ensemble forecast. RASM predicted more sea ice, with an extent further south, in 2020 than in 2019 and this was observed by NSIDC, as well.

PLOS ONE
On the variability of the Bering Sea Cold Pool and implications for the biophysical environment   The latest forecast is presented in Fig 14 for April 15, 2021 and shows an extent that is somewhat similar to 2020. However, the satellite observations show ice extending further south over deep water between 175˚W-180˚. It is important to note that these figures show daily mean sea ice extent and are, therefore, highly influenced by short term atmospheric variability.
Figs 15-17 illustrate the 6-month ensemble forecasts of sea ice extent in Bering Sea using the RASM probabilistic prediction capability. Each ensemble was initialized on February 1st, 2019, 2020, and 2021, respectively. In 2019 and over the Bering Sea only, all the satellite datasets (Bootstrap, NASA Team, and OSI SAF products) show lower ice cover but increased variability of sea ice extent (SIE) relative to the 2009-2018 satellite climatology data, likely due to variable wind forcing. SIE in the Bering Sea decreased through February and reached a winter minimum in March (0.14 ± 0.034×10 6 km 2 ; Fig 15A) before recovering back up to a maximum in late March of 2021. It was one of the lowest sea ice extents in the Bering Sea, which is captured well within the SIE spread of the 31-member ensemble, with the ensemble mean averaging out the temporal variability. The RASM SIE hindcast represents an underestimation, which suggests reduced accuracy of the prescribed atmospheric reanalysis, including near-surface winds. The RASM forecast skill, as measured by the SIE correlation coefficient (CC), was in the range between 0.55 and 0.75 for most ensemble members (74%) and the ensemble mean with their respective normalized standard deviation (SD) values close to those of the satellite SIE SD (Fig 15B). Compared to the satellite climatology, the skill of the ensemble mean forecast represents a gain. Based on these metrics, the RASM hindcast CC was higher (0.89) and its SD was lower (0.51) suggesting a very good agreement in terms of temporal variability but at a reduced magnitude.
In 2020 and 2021, SIE in the Bering Sea increased in February, reached the winter maximum in early March, and then gradually decreased in the following months, with the magnitudes comparable to the 2009-2018 SIE climatology (Figs 16A and 17A). Compared to observations, both the RASM ensemble mean forecast and the RASM hindcast captured the temporal variability but their SIEs were biased low, with only a few ensemble members showing comparable magnitudes. In terms of model forecast skill (Figs 16B and 17B), the ensemble mean CCs were significantly higher than in 2019, at 0.91 and 0.95, and their SDs closer to those in observations, at 0.84 and 1.0, in 2020 and 2021, respectively. While biased low in SIE these skills still represent a clear gain compared to the climatology. Overall, these analyses suggest a potential benefit of the RASM probabilistic forecast capability at time scales of up to 6 months, however it is uncertain how skillful these forecasts would be at longer time scales. Further improvements related to the magnitude of near surface wind speeds and their momentum coupling with the sea ice in the Bering Sea are the subjects of ongoing sensitivity studies.
We show in Fig 18 the predicted bottom water temperature for July 15, 2019, 2020 and 2021 from the same RASM forecasts described above. The predicted CP area of 2019 is the least of the 3 years with a value of 203,000 km 2 , which is only 57% of the long-term mean CP , and the unbiased root-mean-square difference (uRMSE) for the timeseries. The Taylor diagram illustrates a relative model performance in reproducing satellite-derived Bering Sea daily SIE variability, as well as the individual data sets, relative to the average of Bootstrap, NASA NT, and OSISAF products, represented by the black square as the reference where the normalized standard deviation and the correlation coefficient are 1.00 and the unbiased root mean difference is 0.00.

PLOS ONE
On the variability of the Bering Sea Cold Pool and implications for the biophysical environment area during July. This predicted CP area for July 2019 is larger than the July 2018 minimum from the hindcast simulation (113,000 km 2 ). The CP of 2020 is the most extensive of the 3 years with an area of 368,000 km 2 within the 2˚C contour representing the predicted extent of the CP on July 15, 2020. This value is close to the long-term July mean of 359,000 km 2 . The July 15, 2021 CP has an area of 265,000 km 2 in the RASM forecast, which is~26% less than the long-term mean, but still over twice as large as the 2018 value. In a recent study by Kearney et al. [45] it was found that dynamic forecasting had some skill at predicting summer bottom temperatures across the eastern Bering Sea shelf with lead times of up to 4 months, however the skill level was lower when the forecast model was initialized during the fall or early winter. We conclude that while such predictions are now available, it is critical to collect measurements for their ongoing evaluation and potential tuning.